Isochronen in Dresden#
Für die Erklärung bitte hier gucken:
Scenario Beschreibungen:#
Für eine Detailierte Beschreibung bitte in beschreibung-szenarien gucken.
tl;dr
Small: Nur Sidewalks
Medium: Alles außer Große Straßen
Medium-restricted: Sidewalks und Radwege
Large: Nur Straßen keine Fuß und Radwege
Optimal: Alles
cfg = ConfigObj("notebook.ini")
place = "Dresden, Sachsen"
trip_times = [5, 10, 15, 20, 25, 30] # in minutes
travel_speed = 3.6 # in km/h
ox.settings.useful_tags_way += cfr.get_relevant_tags(cfg)
# ox.config(log_console=True, use_cache=True)
Alle Straßen#
Download and prep the street network#
# download the street network
# G_84 = r4r.get_graph(place)
latitude = 51.047993772661066
longitude = 13.74043214154367
G_84 = ox.graph_from_point((latitude, longitude),
dist=1500, network_type='all')
# find the centermost node and then project the graph to UTM
# gdf_nodes = ox.graph_to_gdfs(G_84, edges=False)
# x, y = gdf_nodes['geometry'].unary_union.centroid.xy
# center_node = ox.nearest_nodes(G_84, x[0], y[0])
# center_node = ox.nearest_nodes(G_84, 13.3447448, 50.9162500)
center_node = ox.nearest_nodes(
G_84, longitude, latitude) # Dresden, Rathaus
# center_node = 3930896811
G = ox.project_graph(G_84)
center_node
403203636
# add an edge attribute for time in minutes required to traverse each edge
meters_per_minute = travel_speed * 1000 / 60 # km per hour to m per minute
for u, v, k, data in G.edges(data=True, keys=True):
data['time'] = data['length'] / meters_per_minute
# get one color for each isochrone
iso_colors = ox.plot.get_colors(
n=len(trip_times), cmap='plasma', start=0)
# color the nodes according to isochrone then plot the street network
node_colors = {}
for trip_time, color in zip(sorted(trip_times, reverse=True), iso_colors):
subgraph = nx.ego_graph(G, center_node, radius=trip_time, distance='time')
for node in subgraph.nodes():
node_colors[node] = color
nc = [node_colors[node] if node in node_colors else 'none' for node in G.nodes()]
ns = [15 if node in node_colors else 0 for node in G.nodes()]
fig, ax = ox.plot_graph(G, node_color=nc, node_size=ns, node_alpha=0.8, node_zorder=2,
bgcolor='k', edge_linewidth=0.2, edge_color='#999999', figsize=(16, 16))
Plot the time-distances as isochrones#
How far can you walk in 5, 10, 15, 20, and 25 minutes from the origin node? We’ll use a convex hull, which isn’t perfectly accurate. A concave hull would be better, but shapely doesn’t offer that.
# make the isochrone polygons
isochrone_polys = []
for trip_time in sorted(trip_times, reverse=True):
subgraph = nx.ego_graph(G, center_node, radius=trip_time, distance='time')
node_points = [Point((data['x'], data['y']))
for node, data in subgraph.nodes(data=True)]
bounding_poly = gpd.GeoSeries(node_points).unary_union.convex_hull
isochrone_polys.append(bounding_poly)
# First, plot the network
fig, ax = ox.plot_graph(G, show=False, close=False,
edge_color='#999999', edge_alpha=0.2, node_size=0, bgcolor='k')
# Now, plot each isochrone polygon
for polygon, color in zip(isochrone_polys, iso_colors):
# Extract the exterior coordinates of the polygon
x, y = polygon.exterior.xy
ax.fill(x, y, fc=color, ec='none', alpha=0.6, zorder=-1)
# If there are interiors (holes) in the polygon, plot them too
for interior in polygon.interiors:
x, y = interior.xy
# Fill the hole with the background color (in this case, black)
ax.fill(x, y, fc='k', ec='none', alpha=0.6)
plt.show()
Or, plot isochrones as buffers to get more faithful isochrones than convex hulls can offer#
def make_iso_polys(G, edge_buff=25, node_buff=50, infill=False):
isochrone_polys = []
for trip_time in sorted(trip_times, reverse=True):
subgraph = nx.ego_graph(
G, center_node, radius=trip_time, distance='time')
node_points = [Point((data['x'], data['y']))
for node, data in subgraph.nodes(data=True)]
nodes_gdf = gpd.GeoDataFrame(
{'id': subgraph.nodes()}, geometry=node_points)
nodes_gdf = nodes_gdf.set_index('id')
edge_lines = []
for n_fr, n_to in subgraph.edges():
f = nodes_gdf.loc[n_fr].geometry
t = nodes_gdf.loc[n_to].geometry
edge_data = G.get_edge_data(n_fr, n_to).get(0, None)
if edge_data is None:
continue
edge_lookup = edge_data.get(
'geometry', LineString([f, t]))
edge_lines.append(edge_lookup)
n = nodes_gdf.buffer(node_buff).geometry
e = gpd.GeoSeries(edge_lines).buffer(edge_buff).geometry
all_gs = list(n) + list(e)
new_iso = gpd.GeoSeries(all_gs).unary_union
# try to fill in surrounded areas so shapes will appear solid and blocks without white space inside them
if infill:
new_iso = Polygon(new_iso.exterior)
isochrone_polys.append(new_iso)
return isochrone_polys
# Plotting section
isochrone_polys = make_iso_polys(G, edge_buff=25, node_buff=0, infill=True)
fig, ax = ox.plot_graph(G, show=False, close=False,
edge_color='#999999', edge_alpha=0.2, node_size=0, bgcolor='k')
# Plot each isochrone polygon
for polygon, color in zip(isochrone_polys, iso_colors):
# Extract the exterior coordinates of the polygon
x, y = polygon.exterior.xy
ax.fill(x, y, fc=color, ec='none', alpha=0.6, zorder=-1)
# If there are interiors (holes) in the polygon, plot them too
for interior in polygon.interiors:
x, y = interior.xy
# Fill the hole with the background color (in this case, black)
ax.fill(x, y, fc='k', ec='none', alpha=0.6)
plt.show()
def make_iso_polys_time(G, edge_buff=25, node_buff=50, infill=False, time_window=15):
# Assume trip_times is a list of time values (in minutes) you're interested in
trip_times = [time_window]
isochrone_polys = []
for trip_time in trip_times:
subgraph = nx.ego_graph(
G, center_node, radius=trip_time, distance='time')
node_points = [Point((data['x'], data['y']))
for node, data in subgraph.nodes(data=True)]
nodes_gdf = gpd.GeoDataFrame(
{'id': subgraph.nodes()}, geometry=node_points)
nodes_gdf = nodes_gdf.set_index('id')
edge_lines = []
for n_fr, n_to in subgraph.edges():
f = nodes_gdf.loc[n_fr].geometry
t = nodes_gdf.loc[n_to].geometry
edge_lookup = G.get_edge_data(n_fr, n_to).get(0, None)
if edge_lookup is None:
continue
edge_lookup = edge_lookup.get(
'geometry', LineString([f, t]))
edge_lines.append(edge_lookup)
n = nodes_gdf.buffer(node_buff).geometry
e = gpd.GeoSeries(edge_lines).buffer(edge_buff).geometry
all_gs = list(n) + list(e)
new_iso = gpd.GeoSeries(all_gs).unary_union
if infill:
new_iso = Polygon(new_iso.exterior)
isochrone_polys.append(new_iso)
return isochrone_polys
def convert_coordinates(geom):
"""Convert UTM Polygon/MultiPolygon to Lat/Lon Polygon/MultiPolygon."""
# Define the UTM zone projection and the lat/long projection
utm_zone = "EPSG:32633" # Example for UTM zone 33N, change as per your UTM zone
latlong = "EPSG:4326"
# Change zone accordingly
in_proj = Proj(proj='utm', zone=33, datum='WGS84')
out_proj = Proj(proj='latlong', datum='WGS84')
if geom.geom_type == 'Polygon':
exterior_coords = [transform(in_proj, out_proj, x, y)
for x, y in geom.exterior.coords]
interior_coords = [[transform(in_proj, out_proj, x, y)
for x, y in interior.coords] for interior in geom.interiors]
return Polygon(exterior_coords, holes=interior_coords)
elif geom.geom_type == 'MultiPolygon':
polygons = []
for polygon in geom:
polygons.append(convert_coordinates(polygon))
return MultiPolygon(polygons)
else:
raise ValueError(f"Unsupported geometry type: {geom.geom_type}")
converted_isochrone_polys = [convert_coordinates(
polygon) for polygon in isochrone_polys]
# m = ox.graph_to_gdfs(G_84, nodes=False).explore()
# m
# m = folium.Map(location=[converted_isochrone_polys[0].centroid.y,
# converted_isochrone_polys[0].centroid.x], zoom_start=14)
m = ox.graph_to_gdfs(G_84, nodes=False).explore()
center_point = [latitude, longitude]
# Add each isochrone polygon to the map
for polygon, color in zip(converted_isochrone_polys, iso_colors):
folium.Polygon(
locations=[(lat, lon) for lon, lat in polygon.exterior.coords],
fill=False,
# fill_opacity=0.05,
# fill_color=color,
color=color,
weight=5
).add_to(m)
folium.Circle(
location=center_point,
radius=900, # radius in meters
color='green',
fill=False,
# fill_color='blue',
# fill_opacity=0.05
).add_to(m)
# Add a circle with 2km radius
folium.Circle(
location=center_point,
radius=2000, # radius in meters
color='blue',
fill=False,
# fill_color='blue',
# fill_opacity=0.05
).add_to(m)
# Add a circle with 2km radius
folium.Circle(
location=center_point,
radius=2600, # radius in meters
color='red',
fill=False,
# fill_color='blue',
# fill_opacity=0.05
).add_to(m)
# Display the map
m
Make this Notebook Trusted to load map: File -> Trust Notebook